// cpid.C -- Extracting the calibration constants for the Brho in RIBLL2
// 30Dec, Ni-58 setting, Brho = 4.9896745, brho0*CC=1.6058799
static const int Z = 33;
void cpid(){
  TFile *f = new TFile(getenv("TMP1"));
  TTree *etf = (TTree *)f->Get("etf");
  vector<TH1 *> v;
  TH1 *h; TH2F *h2; TProfile *pfx; TF1 *gaus, *pol1;
  TCut cut = Form("z0>%d.5 && z0<%d.5", Z-1, Z);

  TCanvas *c = new TCanvas("c", "Extracting Brho Correction Constants", 1680, 1050);
  c->Divide(2,2);
  c->cd(1); etf->Draw("bg0>>h1(500, 0.6, 1.)", cut);
  f->GetObject("h1", h); h->Fit("gaus"); gaus = (TF1 *)h->GetFunction("gaus");
  const double bg0m = gaus->GetParameter("Mean");
  c->cd(2); etf->Draw(Form("bg0/%.6f-1:pdcTaUWire[1][0]>>h2(36, -45., 45.,"
  " 500, -0.04, 0.03)", bg0m), cut && "x_f1>=25.&&x_f1<=27.", "colz");
  f->GetObject("h2", h2); h2->ProfileX("pfx_1");
  f->GetObject("pfx_1", pfx); pfx->Fit("pol1", "", "same");
  const double a = pfx->GetFunction("pol1")->GetParameter("p1");
  c->cd(3); etf->Draw(Form("bg0/%.6f-1:x_f1>>h3(52, -52., 52., 500, -0.04, 0.03)",
    bg0m), cut && "pdcTaUWire[1][0]>-2.&&pdcTaUWire[1][0]<2.", "colz");
  f->GetObject("h3", h2); h2->ProfileX("pfx_1");
  f->GetObject("pfx_1", pfx); pfx->Fit("pol1", "", "same");
  const double b = pfx->GetFunction("pol1")->GetParameter("p1");
  c->cd(4); etf->Draw(Form("z0:1.606/bg0*(1+%.8f*pdcTaUWire[1][0]-%.8f*x_f1)*1.04-0.13"
    ">>(1000, 1.9, 2.6, 1000, 10., 40.)", a, b),
    "z0>0.&&x_f1!=-999.&&pdcTaUWire[1][0]!=-999.", "colz");

  c->Print("cpid.pdf");
} // end main function
